Load libraries
library(STutility)
library(ggplot2)
library(ggpubr)
library(magrittr)
library(dplyr)
Assemble spaceranger output files and merge curated meta data
# Mouse brain
samples <- Sys.glob(paths = "../../spaceranger_output/colon/*/filtered_feature_bc_matrix.h5")
imgs <- Sys.glob(paths = "../../spaceranger_output/colon/*/spatial/tissue_hires_image.png")
spotfiles <- Sys.glob(paths = "../../spaceranger_output/colon/*/spatial/tissue_positions_list.csv")
json <- Sys.glob(paths = "../../spaceranger_output/colon/*/spatial/scalefactors_json.json")
infoTable <- data.frame(samples, imgs, spotfiles, json)
infoTable <- cbind(infoTable, arrayID = do.call(rbind, strsplit(infoTable$samples, "/"))[, 5])
curated_metadata <- openxlsx::read.xlsx("../../sheets/curated_RNA_rescue_sample_metadata.xlsx", sheet = 3)
curated_metadata <- setNames(curated_metadata, nm = c("storage_time", "seq_date", "include", "RNA_rescue",
"project", "experimenter", "processer", "comments",
"ID", "paper_id", "RIN", "DV200", "protocol", "source", "arrayID", "spots_under_tissue",
"genes_detected", "fraction_spots_under_tissue",
"median_genes_per_spot", "median_UMIs_per_spot",
"saturation", "reads_mapped_to_probe_set",
"reads_mapped_confidently_to_probe_set",
"reads_mapped_confidently_to_filtered_probe_set",
"reads_mapped_to_genome",
"reads_mapped_confidently_to_genome",
"number_of_panel_genes"))
infoTable <- merge(infoTable, curated_metadata, by = "arrayID")
Load data into a Seurat object
CLN <- InputFromTable(infoTable[c(5:6, 2, 4, 1, 3), ])
## Using spotfiles to remove spots outside of tissue
## Loading ../../spaceranger_output/colon/V11M22-349_A1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../../spaceranger_output/colon/V11M22-349_B1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../../spaceranger_output/colon/V11A20-396_A1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../../spaceranger_output/colon/V11B18-363_C1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../../spaceranger_output/colon/V10S29-108_B1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
## Loading ../../spaceranger_output/colon/V11A20-396_C1/filtered_feature_bc_matrix.h5 count matrix from a 'Visium' experiment
##
## ------------- Filtering (not including images based filtering) --------------
## Spots removed: 0
## Genes removed: 15588
## Saving capture area ranges to Staffli object
## After filtering the dimensions of the experiment is: [21013 genes, 14393 spots]
Check spatial distribution of unique genes
CLN$protocol_array <- paste0(CLN$protocol, " : ", CLN$arrayID)
ST.FeaturePlot(CLN, features = "nFeature_RNA", ncol = 3, label.by = "protocol_array", show.sb = FALSE, pt.size = 1.5)
## Loading required namespace: viridis
Load images
CLN <- LoadImages(CLN, time.resolve = FALSE, xdim = 1e3)
## Loading images for 6 samples:
## Reading ../../spaceranger_output/colon/V11M22-349_A1/spatial/tissue_hires_image.png for sample 1 ...
## Scaling down sample 1 image from 2000x1937 pixels to 1000x968 pixels
## Reading ../../spaceranger_output/colon/V11M22-349_B1/spatial/tissue_hires_image.png for sample 1 ...
## Scaling down sample 2 image from 2000x1937 pixels to 1000x968 pixels
## Reading ../../spaceranger_output/colon/V11A20-396_A1/spatial/tissue_hires_image.png for sample 1 ...
## Scaling down sample 3 image from 1979x2000 pixels to 1000x1011 pixels
## Reading ../../spaceranger_output/colon/V11B18-363_C1/spatial/tissue_hires_image.png for sample 1 ...
## Scaling down sample 4 image from 1979x2000 pixels to 1000x1011 pixels
## Reading ../../spaceranger_output/colon/V10S29-108_B1/spatial/tissue_hires_image.png for sample 1 ...
## Scaling down sample 5 image from 1979x2000 pixels to 1000x1011 pixels
## Reading ../../spaceranger_output/colon/V11A20-396_C1/spatial/tissue_hires_image.png for sample 1 ...
## Scaling down sample 6 image from 1979x2000 pixels to 1000x1011 pixels
Apply rigid transformations to make a rough alignment of the tissue sections.
# Warp transform
CLN <- WarpImages(CLN, verbose = TRUE, transforms = list("1" = list(angle = -90),
"2" = list(angle = -90, mirror.x = TRUE),
"4" = list(angle = -90),
"5" = list(angle = 180),
"6" = list(angle = 180)))
## Creating dummy masks ...Loading masked image for sample 1 ...
## Warping pixel coordinates for 1 ...
## Warping image for 1 ...
## Warping image mask for 1 ...
## Finished alignment for sample1
##
## Loading masked image for sample 2 ...
## Warping pixel coordinates for 2 ...
## Warping image for 2 ...
## Warping image mask for 2 ...
## Finished alignment for sample2
##
## Loading masked image for sample 4 ...
## Warping pixel coordinates for 4 ...
## Warping image for 4 ...
## Warping image mask for 4 ...
## Finished alignment for sample4
##
## Loading masked image for sample 5 ...
## Warping pixel coordinates for 5 ...
## Warping image for 5 ...
## Warping image mask for 5 ...
## Finished alignment for sample5
##
## Loading masked image for sample 6 ...
## Warping pixel coordinates for 6 ...
## Warping image for 6 ...
## Warping image mask for 6 ...
## Finished alignment for sample6
Here, the spots were labeled into three major categories: “mucosa”, “submucosa” and “muscularis”.
The annotations are provided in an Rds file (see below)
CLN <- ManualAnnotation(CLN)
# Add labels
pre_ann <- readRDS("data/CLN_metadata_selections")
CLN@meta.data[rownames(pre_ann), "labels"] <- pre_ann[, "labels"]
Add a new metadata column with a combined protocol and lung ID label
CLN$cln_id <- CLN[[]] %>%
mutate(cln_id = ifelse(paper_id == "1.0", "Sample 1", "Sample 2")) %>%
pull(cln_id)
CLN$protocol_sample <- gsub(pattern = " ", replacement = "_", paste0(CLN$protocol, "_ID: ", CLN$cln_id))
ST.FeaturePlot(CLN, features = "nFeature_RNA", ncol = 3, label.by = "protocol_sample")
Filter out background
CLN <- SubsetSTData(CLN, expression = labels %in% c("mucosa", "submucosa", "muscularis"))
infoTable[c(5:6, 2, 4, 1, 3), 1]
## [1] "V11M22-349_A1" "V11M22-349_B1" "V11A20-396_A1" "V11B18-363_C1"
## [5] "V10S29-108_B1" "V11A20-396_C1"
limits_list <- list("2" = c(x_start = 260, y_start = 500, x_end = 1600, y_end = 1800),
"3" = c(x_start = 170, y_start = 340, x_end = 1400, y_end = 1500),
"4" = c(x_start = 130, y_start = 240, x_end = 1780, y_end = 1560),
"6" = c(x_start = 300, y_start = 280, x_end = 1820, y_end = 1700))
ann_plots <- lapply(c("2", "3", "4", "6"), function(i) {
gg <- cbind(CLN[[]], GetStaffli(CLN)@meta.data)
gg <- subset(gg, sample == i)
dims <- GetStaffli(CLN)@dims[[i]]
p <- ggplot(gg, aes(warped_x, dims$height - warped_y, color = labels)) +
geom_point(size = 1) +
theme_void() +
scale_color_manual(values = c("mucosa" = "#AA4499", "submucosa" = "#DDCC77", "muscularis" = "#CC6677")) +
scale_x_continuous(expand = c(0, 0), limits = c(limits_list[[i]]["x_start"],
limits_list[[i]]["x_end"])) +
scale_y_continuous(expand = c(0, 0), limits = c(dims$height - limits_list[[i]]["y_end"],
dims$height - limits_list[[i]]["y_start"])) +
theme(legend.position = "none")
p <- ggrastr::rasterize(p, layers = "Point", dpi = 300)
return(p)
})
# HE images
plots_HE <- lapply(c("2", "3", "4", "6"), function(i) {
im <- GetStaffli(CLN)@rasterlists$processed[[paste0(i)]]
dims <- GetStaffli(CLN)@dims[[i]]
im <- im[limits_list[[i]]["y_start"]:limits_list[[i]]["y_end"]*0.5, limits_list[[i]]["x_start"]:limits_list[[i]]["x_end"]*0.5]
g <- grid::rasterGrob(im, width = unit(1, "npc"), height = unit(1, "npc"), interpolate = TRUE)
p <- ggplot() +
annotation_custom(g, -Inf, Inf, -Inf, Inf) +
theme_void() +
theme(plot.margin = margin(t = 0, r = 0, b = 0, l = 0))
return(p)
})
p1 <- plots_HE[[1]] +
ann_plots[[1]]
p1
## Warning: Removed 110 rows containing missing values (geom_point).
p2 <- plots_HE[[2]] +
ann_plots[[2]]
p2
## Warning: Removed 161 rows containing missing values (geom_point).
p3 <- plots_HE[[3]] +
ann_plots[[3]]
p3
## Warning: Removed 21 rows containing missing values (geom_point).